Subroutine LUtest

	IMPLICIT NONE
	INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,100)
	INTEGER, PARAMETER :: nspc=3	
	REAL :: fj(nspc),matj(nspc,nspc),inv(nspc,nspc),df(nspc)
	INTEGER :: i,j
	
	open(unit=10,file='ludecout.dat')

	matj=0.0

!	do i=1,nspc
!	  do j=1,nspc
!	  matj(i,j)=i
!	  enddo
!	enddo

	matj(1,1)=1.0
	matj(1,2)=1.0
	matj(1,3)=1.0

	matj(2,1)=0.0
	matj(2,2)=-0.5
	matj(2,3)=1.0

	matj(3,1)=1.0
	matj(3,2)=1.43
	matj(3,3)=0.24

	fj(1)=1.0
	fj(2)=0.5
	fj(3)=1.2

	do i=1,nspc
	write(10,200) (matj(i,j),j=1,nspc)
	enddo

	CALL gauss1(matj,nspc,df,-fj)

	write(10,*) 'inverted matrix:'
	
	do i=1,nspc
	write(10,200) (matj(i,j),j=1,nspc)
	enddo

	write(10,*) 'df:'
	write(10,200) (df(j),j=1,nspc)
!	CALL lut(matj,df,-fj)

!	write(10,*) 'inverse matrix:'
!	do i=1,nspc
!	write(10,200) (matj(i,j),j=1,nspc)
!	enddo


  200	format(3(1pe9.2,2x))
End Subroutine

Subroutine lut(a,x,b)
! LU DECOMPOSITION MADE BY J.C. KIM :
	INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15,100)
	INTEGER, PARAMETER :: nspc = 8

	REAL :: a(nspc,nspc),x(nspc),b(nspc)
	REAL :: factor,sum
	integer :: n,i,j,k
	n=nspc

	do k=1,n-1
	  do i=k+1,n
	    factor=a(i,k)/a(k,k)
	    a(i,k)=factor
	    do j=k+1,n
	      a(i,j)=a(i,j) - factor*a(k,j)
	    enddo
	  enddo
	enddo
	
	do i=1,nspc
	write(*,*) (a(i,j),j=1,nspc)
	enddo
	do i=2,n
	  sum=b(i)
	  do j=1,i-1
	    sum=sum - a(i,j)*b(j)
	  enddo
	  b(i)=sum
	enddo

	x(n)=b(n)/a(n,n)
	do i=n-1,1,-1
	  sum=0.0
	  do j=i+1,n
	    sum=sum+a(i,j)*x(j)
	  enddo
	  x(i)=(b(i)-sum)/a(i,i)
	enddo

	return
end subroutine

Subroutine gauss1(a,n,dx,f)

	IMPLICIT NONE
	INTEGER :: n
	REAL :: a(n,n),dx(n),f(n)

	REAL :: b(n,n),c,d,temp(n),sum
	INTEGER :: i,j,k,m,imax(1),ipvt(n)
	
	b=a
	ipvt=(/(i,i=1,n)/)

	do k=1,n
	  imax=MAXLOC(ABS(b(k:n,k)))
	  m=k-1+imax(1)
	
	  if(m/=k) then
	  ipvt((/m,k/))=ipvt((/k,m/))
	  b((/m,k/),:)=b((/k,m/),:)
	  endif
	  d=1.0/b(k,k)
	
	  temp=b(:,k)
	  do j=1,n
	  c=b(k,j)*d
	  b(:,j)=b(:,j)-temp*c
	  b(k,j)=c
	  enddo
	  b(:,k)=temp*(-d)
	  b(k,k)=d
	enddo

	a(:,ipvt)=b
	
	do k=1,n
	sum=0.0
	  do j=1,n
	  sum=sum+a(k,j)*f(j)
	  enddo
	dx(k)=sum
	enddo
End subroutine
